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ABSTRACT 

Context. Supernova remnants are believed to be a major source of energetic particles (cosmic rays) on the Galactic scale. Since their 
progenitors, namely the most massive stars, are commonly found clustered in OB associations, one has to consider the possibility of 
collective effects in the acceleration process. 

Aims. We investigate the shape of the spectrum of high-energy protons produced inside the superbubbles blown around clusters of 
massive stars. 

Methods. We embed simple semi-analytical models of particle acceleration and transport inside Monte Carlo simulations of OB asso- 
ciations timelines. We consider regular acceleration (Fermi 1 process) at the shock front of supernova remnants, as well as stochastic 
reacceleration (Fermi 2 process) and escape (controlled by magnetic turbulence) occurring between the shocks. In this first attempt, 
we limit ourselves to linear acceleration by strong shocks and neglect proton energy losses. 

Results. We observe that particle spectra, although highly variable, have a distinctive shape because of the competition between ac- 
celeration and escape: they are harder at the lowest energies (index s < 4) and softer at the highest energies (s > 4). The momentum 
at which this spectral break occurs depends on the various bubble parameters, but all their effects can be summarized by a single 
dimensionless parameter, which we evaluate for a selection of massive star regions in the Galaxy and the LMC. 
Conclusions. The behaviour of a superbubble in terms of particle acceleration critically depends on the magnetic turbulence: if B is 
low then the superbubble is simply the host of a collection of individual supernovae shocks, but if B is high enough (and the turbulence 
index is not too high), then the superbubble acts as a global accelerator, producing distinctive spectra, that are potentially very hard 
over a wide range of energies, which has important implications on the high-energy emission from these objects. 

Key words, acceleration of particles - shock waves - turbulence - cosmic rays - supernova remnants 

1 . Introduction account the back-reaction of energetic particles on the shocks). 

However, to ascertain the particle spectrum produced inside 

Superbubbles are hot and tenuous large structures that are the superbubble as a whole, one must also consider important 

formed around OB association s by the powerful winds an d physics occurring between the shocks. Since the bubble interior 

the explosions of massive stars dHigdon & Lingenfeltenl2005l) . is probably magnetized and turbulent, we need to evaluate gains 

They are the major hosts of supernovae in the Galaxy, and and losses caused by the acceleration by waves (a 2nd-order, 

thus major candidates for the production of ene rgetic particles stochastic Fermi process) and escape from the bubble, 

(e.g., lMontmerldfl979l iBvkovlhooU lBuitll2009l and references In this study, we combine the effects of regular acceleration 

therein). Supernovae are indeed believed to be the main con- (occurring quite discreetly, at shock fronts) and stochastic accel- 

tributors of Galactic cosmic rays (along with pulsars and micro- eration and escape (occurring continuously, between shocks), to 

quasars), by means of the diffusive shock acceleration process (a determine the typical spectra that we can expect inside superbub- 

1 st-order, regular Fermi process) occurring at the re mnant's blast bles over the lifetime of an OB cluster. We choose to treat regular 

wave as it goes throug h the interstellar medium dDrurvHl983t acceleration as simply as we can, and concentrate on modeling 



iMalkov & DrurvbOOll) . the relevant scales of stochastic acceleration and escape inside 



Supernovae in superbubbles are correlated in space and time, superbubbles. We present our model in Sect. H give our general 

hence the nee d to in vestigate accel eration by multiple shocks results ln Sect S and P resent s P eclfic ^plications in Sect. H 

dParizotetalJ 12001) . iKlepachet al.l (f2000h developed a semi- Flnall y we dlscuss the limitations of our approach in Sect. |5]and 

analytical model of test-particle acceleration by multiple spheri- P rovlde our conclusions in Sect© 
cal shocks (either wind termination shocks, or supernova shocks 

plus wind external shocks), based on the U miting assumption 2 Model 
of small shocks filling factors. iFerrand et alj d2008l) performed 

direct numerical simulations of repeated acceleration by succes- Our model is based on Monte Carlo simulations of the activity 

sive planar shocks in the non-linear regime (that is, taking into of a cluster of massive stars, in which we embed simple semi- 
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analytical models of (re-)acceleration and escape (described by 
means of their Green functions). To evaluate the average prop- 
erties of a cluster of AT* stars, we perform random samplings 
of the initial mass function (Sect. [XT). For a given cluster, time 
is sampled in intervals d? = 10 000 yr, which is short enough 
to ensure that at most one supernova occurs during that period, 
but by chance for large clusters, and which is long enough to 
consider that regular acceleration at a shock front has shaped 
the spectrum of particles - acceleration is thought to take place 
mostly at early stages of supernova remnant evolution, and in 
a supe rbubble the Sedov phase begins after a few thousands of 
years dParizot et al.ll2004h . Here we do not try to investigate the 
exact extent of the spectrum of accelerated particles: we set the 
lowest momentum (injection momentum) to be p^ a = 10~ 2 m p c 
(which is the typical thermal momentum downstream of a super- 
nova shock) and set the highest momentum (escape momentum) 
to be p max = 10 6 YtipC =i 10 15 eV (which corresponds to the 
"knee" break in the spectrum of cosmic rays as observed on the 
Earth). We note that the theoretical acceleration time from /? m j n 
to /?max (in the linear regime, without escape) is roughly 8 000 yr 
(assuming Bohm diffusion with B = 10^G), which is again con- 
sistent with our choice of df. This corresponds to 8 decades in p, 
at a resolution of a few tens of bins per decade (according to 
Sect.|232J. 

The procedure is then as follows: for each time bin in the 
life of the cluster, either (1) a supernova occurs, and the distri- 
bution of particles evolves according to the diffusive shock ac- 
celeration process, as explained in Sect. 12.21 or (2) no supernova 
occurs, and the distribution evolves taking into account accelera- 
tion and escape controlled by magnetic turbulence, as explained 
in Sect. 12.31 This process is repeated for many random clusters 
of the same size, until some average trend emerges regarding the 
shape of spectra (note that average spectra are not monitored for 
each bin df but in larger steps of 1 Myr). 

In the following, we describe our modeling of massive stars, 
supernovae shocks, and magnetic turbulence. 

2.1. OB clusters: random samplings of supernovae 

We are interested in massive stars that die by core-collapse, pro- 
ducing type lb, Ic or II supernovae, that is of mass greater than 
"Jrnin = 8 m Q , and up to say m max = 120 m Q . These are stars of 
spectral type O (> 20 m ) and include stars of spectral type B 
(4 - 20 m e ). Most massive stars spend all their life within the 
cluster in which they were born, forming OB associations. To 
describe the evolution of such a cluster, one needs to know the 
distribution of star masses and lifetimes. 

The initial mass function (IMF) £ is defined so that the num- 
ber of stars in the mass interval m to m + dm is An = % (m) x dm, 
so that the number of stars of masses between m m j n and m max is 
AT* = J max £ (m) dm . Observations show that £ can be expressed 
as a power law (ISabeterlll955h 

%{m)ccm a , (1) 

with an index of a — 2.30 for massive stars dKroupal2002l) . This 
function is shown in Fig. Q] 

Stars lifetimes can be com puted from stellar evolutio n mod- 
els, and here we use data from lLimongi & Chieffil d2006l) . which 
is plotted in Fig. [2] The more massive they are, the faster stars 
burn their material. A star at the threshold m m ; n = 8 m© has a 
lifetime of fsN.max - 37 Myr, which is also the total lifetime of 
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Fig. 1. Distribution of massive stars masses: the initial mass 
function. For each cluster AT* = 100 stars are randomly chosen 
in the IMF ([]]). The dashed curve represents the experimental 
histogram of masses after A^b = 1000 samples (with resolution 
d log m = 0.05). The dotted curves show 1-, 2-, 3-sigma standard 
deviations over the clusters set. The solid curve is the theoretical 
IMF. 
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Fig. 2. Distribution of massive stars lifetimes (data from 
iLimongi & Chieffil (120061) ). 



the cluster; a star of m max = 120 m G lives only fsN.min - 3 Myr. 
Regarding supernovae, the active lifetime of the cluster is thus 

Ar OB = f SN (m min ) - f SN (ww) - 34 Myr . (2) 



2.2. Supernovae shocks: regular acceleration 
2.2.1 . Green function 

To keep things as simple as possible, we limit ourselves here 
to the test-particle approach (non-linear calculations will be pre- 
sented elsewhere). In the linear regime, we know the Green func- 
tion Gi that links the distribution^] of particles downstream and 
upstream of a single shock according to 

/down(p)= I Gi(p,p ) fu V (po) dpo ; (3) 
Jo 



1 The distribution function f(p) is defined so that the particles num- 
ber density is n = f / (p) Anp 2 dp, where p is the momentum. 
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it reads 



For standard turbulence indices, we obtain 



G\ (p,po) = — — I H(p-po) 
Po \Po) 

where H is the Heaviside function, and 
3r 



si = 



(4) 



(5) 



r- 1 ' 

where r is the compression ratio of the shock. 
2.2.2. Adiabatic decompression 

Around an OB association, particles produced by a supernova 
shock might be reaccelerated by the shocks of subsequent super- 
nova before they escape the superbubble. The effect of repeated 
acceleration is basicall y to harden the spectra (lAchterbergl 1990L 
lMekose&Popelll993l) . 

When dealing with multiple shocks, it is mandatory to ac- 
count for adiabatic decompression between the shocks: the mo- 
menta of energetic particles bound to the fluid will decrease by a 
factor R = r 1 ^ when the fluid density decreases by a factor r. To 
resolve decompression properly, the numerical momentum res- 
olution d log p h as to be significantly smaller than the induced 
momentum shift dFerrand et alj|2008l) . 

2.3. Magnetic turbulence: stochastic acceleration and 
escape 

Particles accelerated by supernova shocks, although energetic, 
might remain for a while inside the superbubble because of mag- 
netic turbulence that scatters them (they perform a random walk 
until they escape). Because of this turbulence, particles will also 
experience stochastic reacceleration during their stay in the bub- 
ble. We present here a deliberately simple model of transport, to 
obtain the relevant functional dependences and order of magni- 
tudes of the diffusion coefficients. The turbulent magnetic field 
8B is represented by its power spectrum W(k), defined so that 

SB 2 oc f m " W(k) dk, where k = 2n/A, A is the turbulence scale, 
and k m i„ (respectively k max ) corresponds to waves interacting 
with the particles of highest (respectively lowest) energy. This 
spectrum is usually taken to be a power law of index q 



W(k) oc k-i , 

normalised by the turbulence level 
(60) 



T] T = 



B 2 + (6B 2 )' 



(6) 



(7) 



2.3.1 . Diffusion scales 

If the turbulence follows Eq. (O, then the space diffusion coeffi- 
cient is given by 



D x (p)=D*x{^-) ' , 
\m p cj 



where we assume that the turbulence spectrum extends suffi- 
ciently for this description to rem ain correct a t the lo west par- 
ticle energies. Using results from ICasse et al.l (120021) obtained 
for isotropic turbulence, one can assume that 



D* oc rT T x B q - 2 A 



9-1 
max 



(9) 



d x (p) 



n T (iO^g) (lOpc) (m' ;C | 1 5 /3 
^(io|g) {wps) [m^) <7 = 3 / 2 



(10) 



Particles diffuse over a typical length scale of x&g = a/6 D x t. 
They are confined within the acceleration region of size jt acc 
as long as x&g (f) < x acc , hence a typical escape time is f esc = 
x acc /6 D x , that is, using Eq. © 



r esc (p) = f* c x | — 
\m p c 



9-2 



where 



oc 



2-q ,1-9 Jl 



Tjj B q A 



"esc ~ " ' L max ■"•ace • 

For standard turbulence indices, we obtain 



■Xp) 



L0 1 



5,o(iOjug) (lOpc) (Wpc) (»ipc) ^ 9 — 5/3 
^(kJ/jg) (lOpc) (iifpc) (h^c) <7 = 3/2 



(11) 



(12) 



■(13) 



Interaction with waves also leads to a diffusion in momen- 
tum. Using results from quasi-linear theory, we can express the 
diffusion coefficient as 

\i 

(14) 



D p (p) = D* p x(m p c) 2 x[-!— 



P 

m„c 



where 



D* oc Tj T B q A 



4-9 l 1 -9 



(15) 



and n is the number density (which determines the Alfven veloc- 
ity together with B). For standard turbulence indices, we obtain 

D„(p) 



g z .cm z .s 



20 ( 10#g) ( lOpc) ( 10- 2 cm--' ) \m p cj 1 ~ 5/3 
m(to^g) (lUpc') (lO- 2 cm- 3 ) ■? = 3/2 



(16) 



2.3.2. Green function 



iBecker et al.l d2006b presented the first analytical expression of 
the Green function G 2 for both stochastic acceleration and es- 
cape that is valid for any turbulence index q e]0, 2[. It is defined 
so that, for impulsive injection of distribution f^t, the distribu- 
tion after time t is 



/end (p 







G 2 (p, po, t) /init (Po) dpo ■ 



Neglecting losses, it can be expressed as 
2-q rp~ Viio? 



(8) G 2 (p,po,i) = 



x exp 



(z + zoMl +0 



(17) 



(18) 



l+q 2 y/zzog) 



2(1-0 l~\2-q' l-t 
zip) = 2 p 2 -" I '{(2-q) ^Dp^ , 

f (/) = exp (2 ( 9 - 2) Dp I ^D*t*^ , 
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Fig. 3. Mean supernovae rate as a function of time. For each clus- 
ter, AT* = 100 stars are randomly chosen in the IMF. The central 
curve represents the experimental mean rate of supernovae after 
A^ob = 1000 samples (with resolution dt = 10 4 years). The top 
curves show 1, 2 and 3 standard deviations over the clusters set. 
The solid curve is the theoretical mean rate of supernovae over 
the cluster's active lifetime ©, i.e. N*/ (fsN.max - ?SN,min)- 



/(At) 




Fig. 4. Distribution of the interval between two successive 
shocks (normalised to the average interval between two super- 
novae). For each cluster, the interval between two successive su- 
pernova is monitored, within the numerical resolution d log At = 
0.05. Colour codes for different numbers of stars A'*, logarith- 
mically sampled between 10 and 500 (purple =10, blue = 27, 
green = 71, orange = 189, red = 500). 



where I (o, x) is the modified Bessel function of the first kind, 
and we recall that D* and ?* sc are defined by Eqs. ( fT31 ) and ( TT2l 
respectively. 

Gi represents the distribution of particles remaining inside 
the bubble. One can also evaluate the rate of particles escaping 
the bubble by dividing G 2 by the escape time given by Eq. ( fTTb : 



(?2,esc (p, P0,t) 



G 2 ip, go, i) _ p 2 q G 2 (p,paJ) 

^esc ip) ^esc 



(19) 



3. Results 

3.1. Distribution of supernovae shocks 

Before presenting the spectra of particles, we briefly discuss the 
temporal distribution of shocks during the life of the cluster, be- 
cause this controls the possibility of repeated acceleration. 



3.1 .1 . Rate of supernovae 

As an illustration of our Monte Carlo procedure, if we count the 
number of supernovae in each time bin [t, t + dt], we can esti- 
mate the mean supernovae rate. The result is sho wn in Fig. [3] In 
agreem ent with the "instantaneous burst" model o flCervino et al.l 
d2000l) . we observe that the distributions of masses and lifetimes 
combine in such a way that, but for a peak at the beginning, the 
rate of supernovae is fairly constant during the cluster's life, and 
can be expressed to a first approximation by 



d«SN N* 



dt 



At, 



N+ x3.10~ 8 yr _1 



(20) 



OB 



where we recall that At^ B is the active lifetime of the cluster, 
given by Eq. ©. 



3.1.2. Typical time between shocks 

Knowledge of the time distribution of supernovae is important to 
acceleration in superbubbles, because, depending on the typical 
interval between shocks, accelerated particles may or may not 
remain within the bubble between two supernovae explosions, 



and thus experience repeated acceleratiorQ. We thus monitor the 
time interval Afj # between two successive supernovae. The re- 
sult is shown in Fig. [4] We note that (1) the most probable time 
interval between two shocks is simply the average time between 
two supernovae A?sn = AfQ B /Af* ; and (2) when time intervals 
are normalised by this quantity, all distributions have the same 
shape independently of the number of stars (apart from very low 
numbers of stars). 

To investigate the probability of acceleration by many suc- 
cessive shocks, we now compute the maximum time Af max that 
a particle has to wait within a sequence of n successive shocks. 
Only particles whose escape time is longer than this value may 
experience acceleration by n shocks. As previously, all distribu- 
tions have the same shape once time intervals are normalised by 
A?sn, and are very peaked, but now the most probable value of 
Af max is a few times longer than the average value (the more suc- 
cessive shocks we consider, the higher the probability of obtain- 
ing an unusually long time interval between any two of them). 
This is summarised in Fig. [3] which shows the most probable 
value of Af max as a function of the number of successive shocks. 
We note that Af max may reach 10 times AfsN> an d that it is an 
imprecise indicator when Af* and n are low. 

3.2. Average cosmic-ray spectra 
3.2.1. General trends 

Proton spectra for clusters of two different sizes inside a typical 
superbubble are shown in Fig. [6] For a given sample, we observe 
a strong intermittency during the cluster lifetime (from blue to 
red), especially at early times. Nevertheless, we clearly see con- 
vergence to an average spectrum as we increase the number N 
of samples (from top to bottom). Comparing results for 10 and 
100 stars (left and right), we see that what actually matters is 
the total number of supernovae N x N+. The limit spectrum ex- 
hibits a distinctive two-part shape, with a transition from a hard 



2 Note that this will also strongly depend on the initial energy of the 
particles: the higher the energy they have gained from one shock, the 
sooner they will escape the bubble, and hence the smaller chance they 
have to be reaccelerated by a subsequent shock. 
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Fig. 5. Maximum time interval between two successive shocks 
in a sequence of n successive shocks (normalised to the aver- 
age interval between two supernovae). Solid curves correspond 
to the most frequent value of Af max (i.e., maxima of the curves 
in Fig. Dotted lines indicate the envelope of the distribution 
(they correspond to a decrease in the maximum value by a fac- 
tor of 10, 100, 1000). Colours code the number of stars N+ in the 
same way as in Fig.|4](note that N+ coincides with the maximum 
number of successive shocks n for which data are available). 

regime (fiat spectrum, of slope s < 4) to a soft regime (steep 
spectrum, of slope s > 4). We also show the escaping spectra in 
Fig. |7] We see that they have the same overall shape, but are a 
bit harder (as highly energetic particles escape first) and of much 
lower normalization. 

Hard spectra at low energies are produced by the combined 
effects of acceleration by supernova shocks (Fermi 1) and reac- 
celeration by turbulence (Fermi 2). Soft spectra at high ener- 
gies are mostly shaped by escape, which preferentially removes 
highly energetic particles. The transition energy is controlled 
by a balance between reacceleration and escape timescales, and 
thus depends on the superbubble parameters. 

3.2.2. Parametric study 

For each cluster, we must define eight parameters N*, r, q, tjt, B, 
^max, n, and x acc , which are more or less constrained. We sam- 
ple the size of the cluster roughly logarithmically between 10 
stars and 500 stars , i.e. N* =10, 30, 70, 200, 500. We consider 
only strong supernova shocks of r — 4. We compare the classi- 
cal turbulence indices q — 5/3 (Kolmogorov cascade, K41) and 
q = 3/2 (Kraichnan cascade, IK65). We consider two different 
scenarios for the magnetic field: if a turbulent d ynamo is oper at- 
ing then B * lOyuG and 6B » B, so that rj T ^ 1 dBvkovl2001l) ; if 
not, then because of the bubble expansion B 1 pG and 5B < B 
(if 5B = B/2, then r\j = 0.2). The external scale of the turbu- 
lence /i max is at least of the order of the distance <i* between two 
stars in the cluster, which, f or a typical OB association radius of 
35 pc (e.g.. lGarmanv|[T994T) . and assuming uniform distribution 
(a quite crude approximation), is 



56 pc 



N 



1/3 



(21) 



which is 26, 12 and 7 pc for 10, 100 and 500 stars respectively. 
However, /l max will be higher if turbulence is driven by super- 
nova remnants, the radius of which increases roughly as 



rsNR - 38 pc 



10 4 yr 



2/5 



in the Sedov-Taylor phase inside a superbubble dParizot et al.l 
120041) . Hence, we consider A nmx = 10, 20,40, 80 pc. We con- 
sider the size of the acceleration region to be of the order of the 
radius of a supernova remnant after our time-step df = 1 000 yr, 
which is x acc = 40 pc according to Eq. d22l . However, in evolved 
superbubbles it might be higher, up to more than 100 pc, so we 
also try 80 pc and 120 pc. The density inside a superbubble is 
always low, and to assess its influence we perform simulations 
with « = 10~ 3 cirT 3 , n — 5 x 10~ 3 cirT 3 , and n = 10~ 2 cm -3 . 
This provides 720 different cases to run. And in each case, we 
have to set the number N of samplings per cluster: convergence 
of average spectra typically requires N* x N - 10 4 , but the gen- 
eral trend is already clear as soon as, N+xN ^ 10 3 , so we simply 
take N — lO^N*. 

We thus had to perform many simulations to explore the pa- 
rameter space. However, interestingly, the effects of the 6 pa- 
rameters relevant to stochastic acceleration and escape q, r\j, B, 
^max, n, and x acc can be summarized b y a single parameter , the 
adimensional number 8* introduced by Becker et all (120061) 



1 



D* t* 

^ p *esc 



which, according to Eqs. dT3T > and dT2b varies as 



For standard turbulence indices, we have 



rfo ( 10//g) ( 1 pc ) ( 40pc ) ( 1 (H cm- 3 ) 9 5/3 

3/2 



n l (i0/*g) (lOpc) (46'pc) ( 



r) q = 



(23) 



(24) 



(25) 



(22) 



For all the possible superbubble parameters considered here, 6* 
ranges from 10~ 4 to 10 +4 . Since we consider only strong super- 
nova shocks of r — 4, the single remaining parameter is the num- 
ber of stars Af* (represented by dots of different colours and sizes 
in subsequent plots), which has a weaker impact on our results. 

To characterize the spectra of accelerated particles, we use 
two indicators, which are plotted in Figs. [8] and [9] We checked 
that the results are independent of the resolution, provided that 
there is at least a few bins per decompression shift. The resid- 
ual variability seen originates mostly in the simulation proce- 
dure itself, which is based on random samplings. In Fig. [8j we 
show the momentum of transition from hard to soft regimes, de- 
fined as the maximum momentum up to which the slope may 
be smaller than a given value (3 or 4 here). Above this momen- 
tum, the slope always remains greater than this value. Below 
this momentum, the slope can be as low as 0, meaning that par- 
ticles pile-up from injection - but we note that it can also hap- 
pen to be > 4 at a particular time in a particular cluster sample, 
since distributions are highly variable. As 9* increases, the tran- 
sition momentum falls exponentially from almost the maximum 
momentum considered (a fraction of PeV) to the injection mo- 
mentum (10 MeV). For rule-of-thumb calculations, one can say 
that the slope can be < 3 up to p = 1/6* GeV. In Fig. [9] we 
show the shallowest slope (corresponding to the hardest spec- 
trum) obtained at a fixed momentum (1 GeV and 1 TeV here). 
As 6* increases, the lowest slope rises from (which is possi- 
ble in the case of stochastic reacceleration) to 4 (the canonical 
value for single regular acceleration in the test particle case). 
As expected, the critical 0* between hard and soft regimes de- 
creases as we increase the reference momentum: the break oc- 
curs around 0* — 10 for p — 1 GeV, and around 8* = 0.01 for 
p = 1 TeV. 
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Fig. 6. Sample results of average spectra of cosmic rays inside the superbubble. The particles spectrum / and its logarithmic slope 
s = d log f/d log p are plotted versus momentum p. The size of the cluster is A 7 * = 10 (left) and A 7 * = 100 (right). The number of 
samplings rises from top to bottom: N = 10, 100, 1000. Other parameters are q — 5/3, B — \QpG, i]j = 1, A m . dx = lOpc, jn; acc = 40pc, 
n = 10 2 cm 3 . 
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Fig. 7. Sample results of average spectra of cosmic rays escaping the superbubble. The particles spectrum / per unit time and its 
logarithmic slope s = d log f/d log p are plotted versus momentum p. The size of the cluster is AT* = 10 (left) and iV* = 100 (right). 
The number of samplings is N = 1000. Other parameters are as in Fig.[6j q = 5/3, B — 10 pG, r\j = 1, A max = 10 pc, x acc = 40 pc, 
n = 10~ 2 cm 3 . 




Fig. 8. Hard-soft transition mo- 
mentum as a function of 6* and 
N*. The transition momentum 
is defined as the momentum up 
to which the particle spectrum 
may have index lower than a 
given threshold: s = 3 at the 
left, s = 4 at the right. 6* is a 
dimensionless parameter defined 
by Eq. l !23t . The number of stars 
=10,30,70,200,500 is coded 
by both dot sizes and dot colours. 
Momentum resolution is 10 bins 
per decompression shift, that is 
50 bins per decade. 



This overall behaviour can be explained by noting that 6* is 
roughly the ratio of the reacceleration time to the escape time. 
Low 6* are obtained when reacceleration is faster than escape, 
allowing Fermi processes to produce hard spectra up to high en- 
ergies, as particles become reaccelerated by shocks and/or tur- 
bulence. In contrast, high 9* are obtained when escape is faster 
than reacceleration, resulting in quite soft in-situ spectra, as par- 
ticles escape immediately after being accelerated by a supernova 
shock. The case 6* — 1 corresponds to a balance between gains 
and losses, in the particular case of which the spectral break oc- 
curs around 10 GeV for s > 4, and around 1 GeV for s > 3. 

4. Application 

4.1. A selection of massive star regions 

We gathered the physical parameters of some well observed 
massive star clusters and their associated superbubbles. The re- 



liability and the completeness of the data were our main selec- 
tion criteria. The parameters useful for our study are: the cluster 
composition (number of massive stars), age, distance, size, and 
the superbubble size and density. We note that we are biased 
towards young objects, since older ones are more difficult to iso- 
late because of their large extensions and sequential formations. 
Information about density is sometimes unavailable. The density 
can span several orders of mag nitude, usually bet ween 10~ 2 and 
10 cm in the central cluster dTorres et al.ll2004j). and betw een 
10 3 and 10 1 cm 3 in the superbubble dParizot et al.ll2004l) . If 
X-ray observations are available, it can be indirectly estimated 
from the thermal X-ray spectrum, given the plasma temperature 
and the column density along the line of sight. In the case of 
a complete lack of data, we accept a mean density of between 
5 x 10~ 3 crrT 3 and 5 x 10~ 2 cirT 3 . Unfortunately, the magnetic 
field parameters can not be directly measured, so that we con- 
sider different limiting scenarios: B — 1 pG and r\j = 0.2 if the 
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Fig. 9. Lowest slope as a function of 
9* and N*. The plots show the lowest 
slope (corresponding to the hardest 
spectrum) reached at a given momen- 
tum: p = 1 GeV at the left, p = 1 TeV at 
the right. 6* is a dimensionless param- 
eter defined by Eq. d23 b . The number 
of stars N* =10,30,70,200,500 is coded 
by both dot sizes and dot colours. 
Momentum resolution is 10 bins per 
decompression shift, that is a 50 bins 
per decade. 



turbulence is low, and B — 10 juG and r\j — 1 if the turbulence is 
high. In each case, we compare our results for turbulence indices 
q = 5/3 and q = 3/2. The maximal scale of the turbulence /l max 
may be taken to be as small as the size of the stellar cluster (espe- 
cially in the case where few supernovae have already occurred), 
or as large as the superbubble itself. 

These quantities are used to estimate the key parameter 
in each of the selected objects using Eq. d25V All the param- 
eters and results are summarised in Table 1 . Before discussing 
the implications of these values, we provide details of the se- 
lected regions in the following two sections, regarding clusters 
found in our Galaxy and in the Large Magellanic Cloud (LMC), 
respectively. 

4.1.1. Galaxy 

We selected 6 objects in the Galaxy. 

- Cygnus region: in this region we identify two distinct 
objects, the clusters Cygnus OBI and OB3, which have 
blown a common superbubble, and the cluster Cygnus OB2. 
We no te that the latter was detected at TeV energies by 
Hegra ( lAharonian et al.ll2005l) as an extended sour ce (TeV 
J2032+4130), and by Milagro dAbdo et all [2007b . as ex- 
tended diffuse emission and at least one source (MGRO 
J2019+37). A supershell was also detected around the 
Cygnus X-ray superbubble, which may have been produced 
by a sequence of starbursts, Cygnus OB2 being the very last. 

- Orion OBI: this a ssociation consists of several subgroups 
dBrown et alJI 19991) . the age of 12 Myrs selected here corre- 
sponds to the oldest one (OB la). 

- Carina nebula: this region is one of the most massive star- 
forming regions in our Galaxy. It contains t wo massive 
stella r clusters, Trumpler 14 and Trumpler 16 dSmith et alj 
120001) . of cumulative size of approximately 10 pc. 

- Westerlund 1: this cluster is very compact although it har- 
bours hundreds of massive stars. The size of the superbubble 
is uncertain, and we assume he re the value of 40 pc reported 
bv lKothes & Dougherty! d2007b for the HI shell surrounding 
the cluster. We no te that Westerlund 1 was detected by HESS 
dOhm et al.ll2009l) . 

- Westerlund 2: the distance to t his cluster remains a m atter 
of debate (see the discussion in lAharonian et al1l2007l) . and 

221, 



we adopt here the estimate of iRauw et al.l d2007D . using it to 



re-evaluate the size obtained by IConti & Crowtherl d2004l) . 
We assume that the giant HII region RCW4 9 of size 100 pc 
is the structure blown by Westerlund 2. iTsujimoto et al.l 
d2007l) provided a spectral fit of the diffuse X-ray emission 
from RCW49, from which we deduce a density ~ 1.5 x 
10~ 3 cm' 3 . We note tha t Westerlund 2 was detected by HESS 
dAharonian et al.ll2007t) . 



4.1.2. Large Magellanic Cloud 

We selected 3 objects in the LMC. All density estimates here 
have been derived from observations of diffuse X-ray emission. 
At the distance of the LMC, these observations usually cover the 
entire structure, so that the density deduced is an average over 
the OB association and the ionised region around it. 

- DEML 192: th is region harbours two massive star clusters, 
LH 51 and 54 dLucke & Hodgelll970b. We deduced the spa - 
tial extensions of both clusters from lOev & Smedlevl d 1998b . 
but these are probably overestimates, because the edges of 
the clusters are not clearly defined. 

- 30 Doradus: this region is quite complex as ca n be seen from 
Chandra observations dTownslev et al.l 120061) . In particular, 
the superbubble extension is difficult to estimate precisely. 
We decid ed to assume th e value given for the 30 Doradus 
nebula bv lWalbornl (1199 ll) . The extension of the star cluster 
may be larger t han the core which harb ours several thou- 
sands of stars (Massey & H unter! 1 1 998b . The core size is 
< 10 pc dMassev & Hunter i l 1998b . it is even estimated to be 
~ 2 pc by Walborn (11991b . The number of massive stars in 
R136 depends on the cluster total mass, estimated to be be- 
tween 5 x 10 4 M Q and 2.5 x 10 5 M Q . Using a Salpeter IMF, 
one finds that 7V*(M > 8M ) =* 400 - 2700. We note that the 
stellar formation in 30 Doradus was sequential and started 
more than 10 Myrs ago dMassev & Hunterlll998b . 

- Nil: this giant HII region harbours several star clusters 
LH9, LH10, LH13, and LH14, probably produced as a se- 
quence of starbursts dWalborn et alj[l999b . Here we mostly 
consider the star cluster LH9 at the center of Nil and the 
shell encompassing it (shell 1 in lMac Low et al.l 1998b . LH10 
is a younger star clu ster with an estimated age of 1 Myr 
dWalborn et al.l 11999b in which no supernova has yet oc- 
curred. The other clusters are less powerful. 
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Table 1. Physical parameters for well observed massive-star forming regions in the Galaxy and in the LMC. 





Cluster 


Superbubble 




Name 




Age 
(iviyrj 


Distance 

(KpCJ 


Size™ 
(pc) 


Size (b) 
(pc) 


Density' ' 
(cm ) 


B=l//G 

q - jf 3 


B=lyUG 

n — ^ n 

q - 3/ L 


B=10//G 

q - 3/3 


B=10//G 

q - 31 L 


Cygnus OB 1/3 


38 (16) 


2-6 (l2) 


1.8 (19) 


24 


80-100 (14) 


0.01? 


5.10 4 
5.10 5 


4.10 2 
3.10 3 


4.10" 
4.10' 


2.10- 2 
1 . 10 ' 


Cygnus OB 2 


750 (5) 


3 _4<i2) 


1.4-1.7' 10 ' 


60 <U) 


450? (5) 


0.02 (5) 


2.10 4 
2.10 5 


9.10' 
7.10 2 


1.10° 
2.10' 


4.10~ 3 
3.10~ 2 


Orion OBI 


30-100 (3) 


12< 2 > 


0.45 (2) 


10 (2) 


140x300 (3) 


0.02-0.03 <4) 


3.10 3 
i 1 n 6 

Z. lu 


4.10' 
7 1 n 3 


3.10"' 
o 1 n 2 

Z. 1U 


1.10 J 

3. IV) 


Carina nebula 


9 


3 (23) 


2.3 (8) 


20 


110 (23) 


0.01? 


2.10 4 

Z. lu 


1.10- 
7 1 n 3 


1.10" 
1 1 n- 


5.10~ 3 

J. lu 


Westerlund 1 


450' 11 


3.3 (1) 


3.9 (13 > 


,0) 


40? <13) 


0.01? 


2.10 3 
^ ID 6 


5.10' 
9 m 4 

Z. 1U 


2.10"' 

J.1U 


2.10- 3 
s in -1 


Westerlund 2 


14 (21) 


2<2D 


g(21, 


1(6) 


100 (21,6) 


0.0015 (24) 


1.10- 
5.10 4 


2.10" 
2.10 2 


9.10" 3 
4.10° 


1.10- 4 

1.10 2 


DEM L192 


135 


3 (20) 


50 


60 (20) 


120xl35 <9) 


0.03 (7) 


3.10 5 
1.10 6 


2.10 3 
4.10 3 


2.10' 
9.10' 


6.10 2 
2.10' 


30 Doradus 


> 400 (22) 


2 (17) 


50 


40 (25) 


200 (25) 


0.09 (27) 


2.10 s 
2.10 6 


1.10 3 

7.10 3 


2.10' 
2.10 2 


6.10 2 
3.10 1 


Nil 


130 


5 (26) 


50 


15x30 (18) 


100xl50 (9) 


0.08 (15) 


9.10 4 
4.10 6 


9.10- 
2.10 4 


8.10" 
4.10 2 


3.10- 2 
8.10-' 



N* is the number of stars with mass > 8M Q (a Salpeter IMF has been assumed, expected for Nil where an index of 2.4 has been used). 
Sizes are either the diameter if the region is spherical, or the large and small semi-axis if the region is ellipsoidal. 
The density is the Hydrogen nuclei density. 

Estimates of 0* are calculated from Eq. d25l l. as explained in Sect. 14. II The range of values of 0* given for each object and for each magnetic 

configuration reflects uncerta i nties in the actual values of bubb le de nsity, accelera t or size and turbulence scal e. 

References : (1) iBrandner et all d2008h. (2) iBrown etalJ d!994h. (3) iBrown et al.1 Jl995h. (4 ) iBurrows et all i 19931). (5) ICash et all dl980l) . 



(6) IConti & Crowthej j2004h. (7) ICooper et al.1 d2004f ). (8) [Davidson & Humphreys! dl997l) . ( 9) iDunne et all d2001h. (10) iHansonl 120031) 



(lDlKnodlsedej (12000). (12 ) iKrTodlseder et alj d200^.. 
d2009l). (16) iMassev etaLl d!995[). (17) iMassev & Huntej d!998^ 
(20) lOev & Smedlevl d!998h 



1998T). 



Kothes & Dougherty! (120071). ( 14) lLozinskava et al. 

(18 ) iNazeetall d2004l). ( 1 9) iNichols-Bohlin & Fesenl 



(21) feauw et all d2007l) . (22) ISelman et all dl999 l). (23) ISmith et al.1 ( 12000) . (24) iTsuiimoto etail 0007 



(25) IWalborrJdl991l) (26) Walborn & Parked d 19921) 7(27) I Wang & Helfandl Jl99ll) 



4.2. Discussion 

In Table 1, we can see that in all cases except for q = 3/2, 
B = 10 j-iG, the critical momentum ~ 1/0+ GeV is in the non- 
relativistic regime. Even if at lower energies the particle distri- 
bution is hard, since pressure is always dominated by relativistic 
particles, one should not expect a strong back-reaction of ac- 
celerated particles over the fluid inside the superbubble, com- 
pared to the case where collective acceleration effects are not 
taken into account. However, if the magnetic field pressure is 
close to equipartitio n with the thermal pressure as suggested by 
iParizot et al.l d2004l) . and provided that the turbulence index q is 
sufficiently low, then the impact of particles on their environ- 
ment has to be investigated. More generally, if q is low enough 
and/or B is high enough, then the superbubble can no longer be 
regarded as a sum of isolated supernovae, but acts as a global ac- 
celerator, producing hard spectra over a wide range of momenta. 

One can wonder how solid these results are, given all the 
uncertainties in the data. In particular, the parameter 0* is very 
sensitive to the accelerator size x max . However x max cannot 
be much lower than a few tens of parsecs (the typical size of 
the OB association) and cannot be much larger than 100 pc 
(the typical size of the superbubble). The maximal scale of the 
turbulence, A max , is even more difficult to estimate, but it also 
ranges between those extrema. Determining precisely these 
spatial scales is complicated by the difficulty of estimating the 
supershell associated with a given cluster, all the more so since 
multiple bursts episodes have occurred (as is likely the case 



in 30 Doradus). In addition, 0* is directly proportional to the 
density, which is not always measured with good accuracy, but 
can usually be rather well constrained to within one order of 
magnitude. The upper and lower values of 6, given in Table 1 
reflect the uncertainties in these three key parameters. In the end, 
we believe that the results presented in Table 1 provide a good 
indication of whether or not collective effects will dominate 
inside the superbubble. Across the range of possible values of 
size and density, the main uncertainty in the critical parameter 
9+ is clearly due to our poor knowledge of the magnetic field 
(how strong the field is, how turbulent it is). It can be seen from 
Table 1 that for a given prescription of the magnetic turbulence, 
the values obtained for both Galactic and LMC clusters are not 
very different from one another. 



5. Limitations and possible extensions 

5. 1 . Regarding shock acceleration physics 

The potentially greatest limitation of our model is its use of a 
linear model for regular acceleration: we have not considered 
the back-reaction of accelerated particles on their accelerator, 
whereas cosmic rays may easily modify the supernova remnant 
shock and therefore the way in which they themselves are accel- 
erated dMalkov & DrunfeOOll) . Since non-linear acceleration is 
a difficult problem, only a few models are a vailable, such as the 
time-asymptotic semi-analytical models of iBerez hko & E llison! 



10 



Ferrand and Marcowith: The spectrum of cosmic rays accelerated inside superbubbles 



d 19991) or iBlasi & Vietril (120051). and the tim e-d ependent nu- 
merica l simulations of Kang & Jones (20071) or iFerrand et al.l 
(120081) . We will include one of these non-linear approaches in 
our Monte Carlo framework in extending our current work. We 
can already note that non-linear effects tend to produce con- 
cave spectra, softer at low energies and harder at high energies 
than the canonical power-law spectrum, and may thus compete 
with reacceleration and escape effects that we have shown to 
have opposite effects. Moreover, non-linearity also occurs re- 
garding the turbulent magnetic field (mandatory for Fermi pro- 
cess to scatter off particles), which remarkably can be produced 
by energetic particles themselves by various instabilities. This 
difficult and still quite poorly under stood process has bee n stud- 
ied by means of MHD simulations ( Jones & Kang 2006), semi- 
analytica l models (lAmato & Blasill2006l) . and Monte Carlo sim- 
ulations dVladimirov et al.ll2006l) . 

Another limitation is that only strong primary supernova 
shocks have been considered (of compression ratio r = 4), but 
since superbubbles are very clumpy and turbulent media, many 
weak secondary shocks are also expected (of r < 4). The com- 
pression ratio r depends on the Mach number Ms according to 



M 2 +: 



where 



M, 



£-50 

cs 



vs 



5000 km/s 



(io 6 k) 



-1/2 



(26) 



(27) 



and us is the shock velocity (of many thousands of km/s in the 
early stages of a remnant evolution) and cj is the speed of sound 
in the unperturbed upstream medium (as high as a few hundreds 
of km/s in a superbubble because of the high temperature T of a 
few millions of Kelvin). In the linear regime, the slope of accel- 
erated particles is determined solely by r according to Eq. ([5]). 
In superbubbles, primary supernova shocks have Ms — 50 and 
already r - 4, leading to s 4; but a secondary shock of say 
Ms — 5 has only r 3, leading to s 4.5. We note that al- 
though weaker shocks produce softer individual spectra, being 
more numerous they may help to produce hard spectra by re- 
peated acceleration, so that their net effect is not obvious. To 
begin their investigation, we added a weak shock at each time- 
step immediately following a supernova (except if another su- 
pernova occurs at that moment), of compression ratio randomly 
chosen between 1.5 and 3.5. For regular acceleration alone, im- 
portant differences are seen between simulations including only 
strong shocks, or only weak shocks, or both. But once combined 
with stochastic acceleration and escape, these differences are no 
longer evident. We have repeated our 720 simulations at medium 
resolution and observed that our two indicators (momentum of 
transition and minimal slope) remain globally unchanged. The 
shape of cosmic-ray spectra thus seems to be mostly determined 
by the interplay between reacceleration and escape, acceleration 
at shock fronts acting mostly as an injector of energetic parti- 
cles. We note that, before supernova explosions, the winds of 
massive stars, not explicitly considered in this study, may also 
act as injectors in the same way, as they have roughly the same 
mechanical power integrated over the star lifetime. 

Finally, one may question our particular choice of stellar evo- 
lution models, but we believe that possible variations in the ex- 
act lifetime of massive stars would bring only higher order cor- 
rections to the general picture that we have obtained. We also 
note that we have implicitly considered that stars are born at the 
same time, and then evolve independently, while in reality star 



formation may occur through successive bursts within a same 
molecular cloud, which could be sequentially triggered by the 
first explosions of supernovae. Another possible amendment to 
our model is that stars of mass greater than 40 solar masses may 
end their life without collapsing, and thus without launching a 
shock. We have repeated our 720 simulations at low resolution 
considering the occurrence of supernovae only for m < 40 m Q , 
and checked that our two indicators remain globally unchanged. 
This seems consistent with the shape of the IMF (there are very 
few stars of very high mass) and the shape of star lifetimes (stars 
of very high mass have roughly the same lifetime). 

5.2. Regarding inter-shock physics 

We use an approximate model of stochastic acceleration, be- 
cause of the use of relativistic formulae a nd the neglec t of en - 
ergy losses, to be able to use results from iBecker et alj d2006l) . 
However, we note that, in terms of stochastic acceleration, the 
relativistic regime is reached when m p v » m p \>A, where v is the 
particle velocity and va the Alfven velocity 



B 



va 



2.10 7 cm.s" 1 



B 



10 //G 



(lO- 2 cm- 3 ) 



(28) 



and in a superbubble this condition is met for p » 1 MeV, since 
Va/c - 10~ 3 . Although we could of course implement more in- 
volved models of transport, we emphasize that our main objec- 
tive was to find the key dependences of the problem, and we 
have shown that it is mainly contr olled by the param eter 0* . 
Regarding losses, the formalism of IBecker et al.l d2006l) allows 
for systematic losses, but for mathematical convenience these 
are supposed to occur at a rate oc p q ~ x , which can describe 
Coulomb losses only in the very special case of q — 2. But pro- 
ton losses above 1 GeV are dominated by nuclear interactions 
(lAharonian & Atovanl 19961 ) with a typical lifetime of 6.10 7 yr/n, 
where n is the density in crrT 3 , which is far longer than the su- 
perbubble lifetime given the low density n < 10 -2 cm' 3 (but 
this might become a concern when cosmic rays reach the parent 
molecular clouds where n > 10 2 cm 4 ). At very low energies 
(around the MeV), ionization losses might also be important and 
compete with stochastic reacceleration. 

Finally, we note that most parameters are time-dependent, 
and might become considerably different at late stages. For com- 
pleteness, we have performed our simulations until the explo- 
sion of the longest lived stars, but over tens of millions of years 
the overall morphology and properties of the superbubble might 
change substantially as it interacts with its environment. As 
long as the typical evolution timescale of relevant parameters 
is longer than our time-step dt = 10 000 yr, their variation can 
be taken into account simply by varying the value of 9* accord- 
ingly. Otherwise, direct time-depende nt numerical simulations 
similar to those o f IFerrand et ail d2008l) will be necessary. 



6. Conclusions 

Our main conclusions are as follows: 

1. Cosmic-ray spectra inside superbubbles are highly variable: 
at a given time they depend on the particular history of a 
given cluster. 

2. Nevertheless, spectra follow a distinctive overall trend, pro- 
duced by a competition between (re-)acceleration by regular 
and stochastic Fermi processes and escape: they are harder at 
lower energies (s < 4) and softer at higher energies {s > 4), 
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shapes that are in agreement with the results of lBvkovl ( l2001l) 
based on different assumption^. 

3. The momentum at which this spectral break occurs critically 
depends on the bubble parameters: it increases when the 
magnetic field value and acceleration region size increase, 
and decreases when the density and the turbulence external 
scale increase, all these effects being summarized by the sin- 
gle dimensionless parameter 8* defined by Eq. (1231 . 

4. For reasonable values of superbubble parameters, very hard 
spectra (s < 3) can be obtained over a wide range of ener- 
gies, provided that superbubbles are highly magnetized and 
turbulent (which is a debated issue). 

These results have important implications for the chemistry 
inside superbubbles and the high-energy emission from these ob- 
jects. For instance, in the superbubble Perseus OB 2 there is ob- 
serva tional evidence of intense spallation activity dKnauth et al.l 
120001) attributed to a high density of low-energy cosmic rays, 
but EGRET has not detected 7r°-decay radiation, which places 
strong limits on the density of high-energy cosmic rays. This 
is consistent with the shape of the spectra obtained in this work. 
We are thus looking forward to seeing how new instruments such 
as Fermi and AGILE will perform on extended sources such as 
massive star forming regions, which have recently been estab- 
lished as very high-energy sources. In that respect, we make a 
final comment that the high intermittency of predicted spectra 
might explain the puzzling fact that some objects are detected 
while others remain unseen. 

Acknowledgements. The authors would like to thank Isabelle Grenier and 
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